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Abstract 

Siphons in a chemical reaction system are subsets of the species 
that have the potential of being absent in a steady state. We present 
a characterization of minimal siphons in terms of primary decompo- 
sition of binomial ideals, we explore the underlying geometry, and we 
demonstrate the effective computation of siphons using computer al- 
gebra software. This leads to a new method for determining whether 
given initial concentrations allow for various boundary steady states. 

Keywords: chemical reaction systems, siphon, steady state, mono- 
mial ideal, binomial ideal, primary decomposition 

1 Introduction 

In systems biology, a population model or chemical reaction system is said 
to be "persistent" if none of its species can become extinct if all species are 
present at the initial time. Those subsets of the species that can be absent 
in steady state are called "siphons." Angeli et al. [5] suggested the concept 
of siphons to study the long-term behavior of dynamical systems that model 
chemical reactions. In terms of the dynamics, a siphon is the index set of a 
forward-invariant face of the positive orthant. Any boundary steady state 
must lie in the interior of such a face. Hence, to investigate the trajectories, 
it is useful to list all minimal siphons. The present paper offers an algebraic 
characterization of siphons, and it shows how this translates into a practical 
tool for computing siphons. 

Following [H [9] , we represent a chemical reaction network as a directed 
graph G whose nodes are labeled by monomials and whose edges correspond 
to reactions. A siphon of G is a non-empty subset Z of the variables such 
that, for every directed edge m — > m' in G, whenever one of the variables in 
the monomial m' lies in Z then so does at least one of the variables in m. 
In Section 2 we relate this definition to the description of siphons given 
in O [7], we review the underlying dynamics, and we discuss its meaning 
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in terms of polyhedral geometry. Our algebraic approach is presented in 
Section 3. Theorem 13.11 expresses the minimal siphons of G in terms of the 
primary decomposition of a binomial ideal associated to G. If the directed 
graph G is strongly connected then the ideal encoding the minimal siphons 
is generated by squarefree monomials. In Theorem l3.2l and Algorithm 13. 61 we 
explain how to compute the relevant (stoichiometrically compatible) siphons 
for any set of initial conditions. In particular, a chemical reaction system 
without relevant siphons has no boundary steady states, and this property 
is sufficient for proving persistence in many systems [3 [20]. In Section 4, 
we demonstrate that the relevant computations can be performed effectively 
using computer algebra software, such as Macaulay 2 [15j. 

In the remainder of the Introduction we present three examples from the 
systems biology literature, with the aim of illustrating our algebraic repre- 
sentation of chemical reaction networks and the computation of siphons. 

Example 1.1. We consider a receptor-ligand dimer model, which is ana- 
lyzed by Chavez in her thesis ^ §7.2] and by Anderson [21 Example 4.1]: 
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Note that the reaction A^C AD is usually denoted hy 2A + C ^ A + D. 
The biochemical species are as follows: the species A denotes a receptor, 
B denotes a "dimer" state of A (two receptors joined together), and C de- 
notes a ligand that can bind either to A (to form D) or to B (to form 
E). There are three minimal siphons, {A,B,E^, {A,G,E}, and {C, 
which correspond to the minimal primes of the monomial ideal of the com- 
plexes (j4^C, AD, E, BC). We will return to this example in Section m □ 

Example 1.2. The following enzymatic mechanism was analyzed by Siegel 
and MacLean |20j, and also by Chavez |6], Example 4.6.1]: 
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The species are S (a substrate), E (an enzyme), P (a product), I (an uncom- 
petitive inhibitor), and intermediate complexes Q and R. Here the graph 
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consists of two strong components, and we encode it in the binomial ideal 
{SE -Q, Q- PE, QI -R) + {EPQRS). The radical of this ideal equals 

{E, Q, R) n (/, R, ES-Q, P-S) n {P, Q, R, S) . 

By Theorem l3.lt the minimal siphons are the variables in these prime ideals. 
Thus the minimal siphons are {E, Q, R}, {I, R}, and {P, Q, R, S}. □ 

Example 1.3. Here is the network for a basic one-step conversion reaction: 

SqE ^ X ^ PE 
PF ^ Y ^ SqF . 

The enzyme E helps convert a substrate into a product P, and a second 
enzyme F reverts the product P back into the original enzyme Sq; these are 
also called "futile cycles" [H |T7] . Such reactions include phosphorylation and 
de-phosphorylation events, and they take place in MAPK cascades. This 
network has three minimal siphons: {E, X}, {F, Y}, and {P, So, X, Y}. 
To see this algebraically, we form the binomial ideal 

{ESo-X, X{EP-X), FP-Y, Y{FSo-Y), EFPSqXY) . 

This ideal corresponds to Tg in Theorem 13. H and it has six minimal primes: 

{E, X, FSo-Y, P-So), {F, Y, P - So, SqE-X), {P, So, X, Y), 
{E, X, Y, F), {E, X, P, Y), and {F, So, X, Y) . 

The three minimal siphons arise from the first three of these six primes. □ 

2 Reaction networks, siphons, and steady states 

A chemical reaction network is defined by a finite labeled directed graph G 
with n vertices. The i-th vertex of G is labeled with a monomial c^* = 
c^'^Cg'^ • • • Cs" in s unknowns ci, . . . ,Cs, and an edge {i,j) is labeled by a 
positive parameter Kij. This graph defines the ordinary differential equations 

ft = w 

where ^(c) = (c^^, c^^, . . . , c^") is the row vector of the monomials, Y = 
{Uij) is the n x s-matrix of exponent vectors of the n monomials, and is 
the n X n-matrix whose off-diagonal entries are the Kij and whose row sums 
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are zero (i.e. minus the Laplacian of G). The equations ([T]) are those of 
mass- action kinetics, although the concept of a siphon is independent of the 
choice of kinetics. In order for each chemical complex c^' to be a reactant or 
product of at least one reaction, we assume that G has no isolated points. 
For a complex c^^ and for a G [s], we write Ca\c^'^ ("Ca divides c^^*") if yia > 0; 
in other words the i-th complex contains species a. If the i-th complex does 
not contain species a, then we write Ca f c^'. 

A non-empty subset Z of the index set [s] := {1, 2, . . . , s} is a siphon if 
for all z Z and all reactions c^' — > c^^ with c^|c^^ , there exists a Z 
such that Ca|c^'. Siphons were called "semilocking sets" in [2l|3]. Note that 
the set of siphons of G does not depend on the choice of parameters Kij. 

With any non-empty subset Z C [s] we associate the prime ideal 



in the polynomial ring Q[ci, C2, . . . , c^]. Recall (e.g. from [8]) that the variety 
of ^z, denoted by V{^z), is the set of points x € such that /(x) = 
for all polynomials / G ^z- Thus, the non- negative variety V>q{^z) is the 
face of the positive orthant M>q defined by all Z-coordinates being zero. 

Proposition 2.1. A non-empty subset Z of [s] is a siphon if and only if 
V>o{^z) is forward-invariant with respect to the dynamical system 

Proof. This is the content of Proposition 2 in Angeli et aZ. jS]. □ 

In Example I l.H the dynamical system ([T]) takes the explicit form 

dA/dt = - {ki2 + 2ku)A^G + {k2i - k23)AD + K32E + Ik^BC 
dB/dt = Ki4,A^G - {k4i -\- Ki3)BG + K34E 
dG/dt = - kuA'^C + K21AD - K43BC + K34E 

dD/dt = KuA^C - {k21 + K23)AD + K32E 

dE/dt = K23AD + K43BG - {k32 + H3i)E . 

This is a dynamical system on M>q. Each of the three minimal siphons 
{A, B, E}, {A, G, E}, and {C, D, E} defines a two-dimensional face of M>q. 
For example, ^>o(*P{a,_b,_e}) is the face in which the coordinates A, B, and 
E are zero and C and D are non-negative. The minimality of the three 
siphons implies that no face of dimension three or four is forward- invariant. 

We next collect some results relating siphons to boundary steady states, 
that is, non-negative steady states of ([TJ having at least one zero-coordinate. 
These connections are behind our interest in computing siphons. See [2l[3l[5] 
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for details on how siphons relate to questions of persistence (the property 
that positive trajectories of ([T]) have no accumulation points on the boundary 
of the orthant K>q). We first show that a boundary steady state necessarily 
lies in the relative interior of a face V>o{^z) indexed by a siphon Z. 

Lemma 2.2. Fix a reaction network G, and let j be a point on the boundary 
of the positive orthant M>q with zero coordinate set Z := {i [s] : 7^ = }. 
If 'J is a boundary steady state of (OP, then the index set Z is a siphon. 

Proof. Assume that c^jc^J for some species z £ Z and some complex c^^ of 
G. Let I index complexes that react to c^^ but do not contain the species z: 

T := {z G [n] : c^* — > c^J is a reaction of G and Cz /fc^'}- 



J^/^iiyi.T^^ = 0, (2) 

where the second equality holds because 7 is a steady state. The summands 
of dl]) are non-negative, so we have 7^* = for all z G X. Thus if i € Z there 
exists a, G [s] with 7^^ = (so, Oj G Z) and hence Cajc^'- □ 

A similar result holds for boundary w-limit points (accumulation points) 
of a trajectory; see [2[ [5] or [31 Theorem 2.13]. We are interested in the 
dynamics arising from some initial condition c^^^ G M?^q, so we restrict our 
attention to polyhedra, called invariant polyhedra, of the following form: 

^c(o) := (c(°)+Lstoi) n RU . (3) 

Here Lstoi := spanjyj — yi : c^* cP^ is a reaction} is the stoichiometric 
subspace in M*. For any c and any k, the right hand side vector ^'(c) • • Y 
of the dynamical system ([1]) lies in Lstoi , and therefore the polyhedron V^{o) 
is forward-invariant with respect to ([T]). For any index set W C [s], let 

Fw ■■= {xeV^m : Xi = OifieW} = F>o(*Ph/) nP,(o) 

denote the corresponding (possibly empty) face of V^(o) ■ All faces of V^m 
have this form; see (S] §2.3] for further details. Lemma 12.21 implies the 
following: Given an invariant polyhedron V^m , if all siphons Z yield empty 
faces, Fz = 0, then V^{o) contains no boundary steady states. In Theorem 
13.51 we shall present an algebraic method for deciding when this happens. 

We now examine the case when the chemical reaction network is strongly 
connected, i.e., between any two complexes there is a sequence of reactions. 



Then we have 
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Lemma 2.3. Assume that G is strongly connected. Then a point 7 € M>g 
is a boundary steady state if and only if Z = {i £ [s] : 'ji = 0} is a siphon. 

Proof. The forward implication is Lemma 12.21 Now let 7 be a boundary 
point whose zero-coordinate set Z is a siphon. Because G is strongly con- 
nected, all complexes c^' evaluated at 7 are zero (7^' = 0), and hence, each 
monomial that appears on the right-hand side of ([T]) vanishes at c = 7. □ 

Prom a polyhedral geometry point of view, Lemma [2.31 states the follow- 
ing: For strongly connected reaction networks G, any face of an invariant 
polyhedron V^^o) either has no steady states in its interior or the entire face 
consists of steady states. We shall see now that a similar result holds for 
toric dynamical systems. Recall (from ^) that ([T]) is a toric dynamical sys- 
tem if the parameters Kij are such that ^'(c) • Af^ = has a positive solution 
c € MS^Q (which is called a complex-balancing steady state). The following 
result concerns the faces of invariant polyhedra of toric dynamical systems. 

Lemma 2.4. Let c^^^ G be a positive initial condition of a toric dy- 

namical system. Then a face Fz of the invariant polyhedron V^^o) contains 
a steady state in its interior if and only if Z is a siphon. 

Proof. This is derived in [21 Lemma 2.8], albeit in different language. □ 

3 Binomial ideals and monomial ideals 

In what follows we characterize the minimal siphons of a chemical reaction 
network G in the language of combinatorial commutative algebra [18] . It will 
be shown that they arise as components in primary decompositions. For any 
initial conditions c^^\ we characterize those siphons that define non-empty 
faces of the invariant polyhedron V^(o) ■ In the next section we shall see that 
these results translate into a practical new method for enumerating siphons. 

Throughout this section we fix the ring R = Q[ci, . . . , Cs]/ {ciC2 . . . Cs) ■ 
This is the ring of polynomial functions with Q-coefficients on the union of 
the coordinate hyperplanes in M^. All our ideals will live in this ring. 

With a given network G we associate the following three ideals in R: 

= ( c^' • (c^J - c^') : c^' c^^ is a reaction of G), 
:Sq = c^j - c^' : c^' c^J is a reaction of G ) , 

tig = {^{c)) = {cy\cy\...,cy"). 

Thus Tg encodes the directed edges, and encodes the underlying undi- 
rected graph. These are pure difference binomial ideals [12l [13], while 9JTg 
is the monomial ideal of the complexes. The following is our main result. 
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Theorem 3.1. The minimal siphons of a chemical reaction network G are 
the inclusion-minimal sets {i G [s\ : c, G ^} where ^ runs over the minimal 
primes of%G- If each connected component of G is strongly connected then 
1g con be replaced in this formula by the ideal ^g ■ Moreover, if G is strongly 
connected then %g can be replaced by the monomial ideal TIq ■ 

Proof. The complex variety Vc(Tg) consists of all points 7 € having at 
least one zero coordinate and satisfying 7^' • (7^^ —jVi'j = Q for all reactions. 
We first claim that our assertion is equivalent to the statement that the 
minimal siphons are the inclusion-minimal sets of the form {i G [s] : 7i = 0} 
where 7 runs over Vc(Tg)- Indeed if ^ is a minimal associated prime of Tg, 
let 7 € {0, 1}* be defined by 7j = 1 if and only if Cj ^ *p. It follows that 
{i e [s] : Q G qj} = {i e [s] : 7, = 0} and 7 G Vc^i) C Vc{1g)- Conversely, 
if 7 G Vc('Jg'), then 7 G Vc(*P) for some minimal associated prime and 
so we have the containment {i G [s] : Cj G C {i G [s] : 7^ = 0}. If, 
furthermore, the set {i G [s] : 7i = 0} is minimal among those defined by 
7' G Vc('Xg'), then by above it must follow that the containment is in fact 
equality. 

Next, if 7 is in Vc(Tg) then we can replace 7 by the 0-1 vector 5 defined 
by <Ji = if 7i = and 5j = 1 if 7^ ^ 0. This non-negative real vector has the 
same support as 7 and lies in the variety of %g- Hence our claim is that the 
minimal siphons are the inclusion-minimal sets of the form {i G [s] : 6i = 0} 
where 6 runs over 1^{o,i}('^g)- But this is obvious because 6^^ ■ {6^^ — 6^^) = 
if and only if S^^ = implies 6^^ =0. 

Now, the minimal associated primes of Tg depend only on the radical of 
Tg, so we can replace Tg by any other ideal that has the same radical. If the 
components of G are strongly connected then the complex c^' can produce 
c^i if and only if c^J can produce and in this case both c^' • (c^j — c^*) 
and c^J • (c^* — c^^ ) are in Tg- Hence the radical of Tg contains the binomial 
(.Vi — cVj ^ and we conclude that Tg and Zg have the same radical. 

Finally, SOTg is a monomial ideal, and associated primes of a monomial 
ideal are of the form '^z for some Z d [s\. It is straightforward to see that 
if G is strongly connected, '^z contains OJTg if and only if Z is a siphon. □ 

When analyzing a concrete chemical reaction network G, one often is 
given an initial vector c^^^ G MS^q for the dynamical system ([T]), or at least 
a subset of MS^q that contains c^^\ A siphon Z C [s] of G is called 
c^^^ -relevant if the face Fz of the invariant polyhedron V^io) is non-empty. 
In other words, if Z is c'-''^ -relevant, then there exists a boundary point 
that is stoichiometrically compatible with c^'^^ and has zero-coordinate set 
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containing Z. For any subset 17 of MS^q, we say that Z is il-relevant if it is c^^^- 
relevant for at least one point c^^^ in ft. Finally we call a siphon relevant if it 
is RS^Q-relevant. Relevant siphons are also called "critical" siphons [5j, and 
non-relevant siphons are also called "stoichiometrically infeasible" siphons 
[3] and "structurally non-emptiable" siphons [S]. We next explain how to 
enlarge the ideals Tg, 3g) and Tic so that their minimal primes encode only 
the siphons that are relevant. 

We recall that the stoichiometric subspace Lstoi of M'' is spanned by all 
vectors Uj — yi where c^' — ?> c^^ is a reaction in G. Its orthogonal complement 
-^cons := (^stoi)^ is the space of conservation relations. Let Q denote the 
image of the non-negative orthant M>q in the quotient space M*/Lgtoi — 
-^cons- Thus Q is a convex polyhedral cone and its interior points are in 
bijection with the invariant polyhedra V^io)- Further, Q is isomorphic to the 
cone spanned by the columns of any matrix A whose rows form a basis for 
-^cons- This isomorphism is given by the map 

: Q ^ l^aiUi : Qi,a2, . . . ,as > ^ , ^ i-^ gjQj , 

where q = {qi,q2, • • • , 9s) G ^>o and ai, a2, • • • , are the columns of the 
matrix A. For simplicity, we identify the cone Q with the image of <j)A- 
A subset F of [s] = {1, 2, . . . , s} is called a facet of Q if the corresponding 
columns of A are precisely the rays lying on a maximal proper face of Q. 
Any maximal proper face of Q also is called a facet. The list of all facets of 
Q can be computed using polyhedral software such as polymake [H]. 

We represent the facets of Q by the following squarefree monomial ideal: 

= Pi (c, :i0F> = fl ^p.. 

F facet of Q F facet of Q 

Each vertex of an invariant polyhedron V^n)) is encoded uniquely by its 
support V , which is a subset of [s]. Consider the squarefree monomial ideal 

^cC) ~ ( W • ^ encodes a vertex of V^m ) . 

The distinct combinatorial types of the polyhedra V^m determine a natural 
chamber decomposition of the cone Q into finitely many smaller cones: if 
two polyhedra V^{o) and V^(o) correspond to points in such a chamber of 
the decomposition, then the polyhedra have the same set of supports V 
of their vertices. For an example see Figure [TJ In the context of chemical 
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reaction networks, the same chamber decomposition appeared in recent work 
of Craciun et al. [10]. Specifically, its chambers were denoted Si in [lOl §2.1]. 

The ideal ^^{o) depends only on the chamber that contains the image of 
c^'^\ For any subset C M?^q, we take the sum of the ideals corresponding 
to all chambers that intersect the image of in Q. That sum is the ideal 

*Bq = ( Y\ Cj : V encodes a vertex of V^io) for some c^^^ GO,), 
iev 

The above ideals are considered either in the polynomial ring Q[ci, . . . ,Cs] 
or in its quotient R = Q[ci, . . . , Cg]/ (ciC2 • • • Cg), depending on the context. 

Let Ti and T2 be two arbitrary ideals in R. Recall (e.g. from [8]) that the 
saturation of Ti with respect to T2 is a new ideal that contains Ti, namely, 

SatCri,T2) = (Ti : T^) = {f £ R : f ■ (T2)™ C Ti for somem G Z>o}. 

Here, we shall be interested in the following nine saturation ideals: 

SatCTc,®), Sat (Tg, 53^(0)), SatCTc^So), 
Sat(aG,«), Sat (^G, 55,(0)), Sat(aG,5Sn), (4) 
Sat(mtG,'B), Sat (OKg, 53,(0)), Sat(aJtG, 5Sn). 

The following theorem is a refinement of our result in Theorem 13.11 

Theorem 3.2. The relevant minimal siphons ofG are the inclusion-minimal 
sets {i € [s] : Ci G *p} where *p runs over minimal primes from The ide- 
als in the first, second, and third columns yield relevant siphons, c^^^ -relevant 
siphons, and Q-relevant siphons, respectively. The ideals in the first row are 
for all networks G, those in the third row for strongly connected G, and those 
in the middle row for G with strongly connected components. 

Proof. The variety of the ideal Sat(Ti,T2) is the union of all irreducible 
components of the variety V{1-i) that do not lie in V{%2)- The result now 
follows from Theorem 13.11 and the following observations. The non-negative 
variety V>o{^) consists of all points in M>q whose image modulo Lgtoi lies in 
the boundary of the cone Q. Thus, for a minimal siphon Z, the image of the 
variety V>q{^z) is in the boundary of Q if and only if Z is not relevant. More 
precisely, the image of V>o{^z) is in the interior of the subcone spanned by 
{oj : i ^ Z}, so there exists a facet of Q that contains the subcone if and 
only if Z is not relevant. Therefore, any irreducible component of V{2g) 
(or V^Zg) or V{d)lG)) defines a non-relevant siphon Z if and only if it lies 
in Vl^ip") for some facet F of Q, which is equivalent to lying in V{^). 
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Next, the variety y>o(5S^(o)) is the union of all faces of the orthant M>q 
that are disjoint from the invariant polyhedron ^^(o) • So, for a minimal 
siphon Z, the ideal does not contain *B^(o) if and only if there exists a 
vertex of V^^o) whose zero-coordinate set contains Z, which is equivalent to 
the condition that the face Fz of the polyhedron is non-empty. Hence, any 
component of V{^g) (o^ V{1g) or V{TIg)) that defines a minimal siphon Z 
lies in V (53^(o) ) if and only if Z is not relevant. Finally, the variety F>o(5Sn) 
is the intersection of the varieties y>o(5S^{o)) as c^^^ runs over 0. □ 

Example 3.3. In Example 1 1 . 2 1 and Example 11.31 Q is the cone over a trian- 
gle, and the three minimal siphons are precisely the facets of that triangular 
cone. Thus, there are no relevant siphons at all. This is seen algebraically 
by verifying the identities Sat{^G,^) = (1) and Sat(TG,*B) = (1). □ 

We now discuss the case when a network has no relevant siphons, by 
making the connection to work of Angeli et al. [5] , which focuses on chemical 
reaction networks whose siphons Z all satisfy the following condition: 

(■^) there exists a non-negative conservation relation I G Lcons H IR>o 
whose support supp(/) = {i G [s] : > 0} is a subset of Z. 

Recall that Angeli et al. call siphons satisfying this property "structurally 
non-emptiable" [5l §8]. Note that the property (-^) needs only to be checked 
for minimal siphons in order for all siphons to satisfy the property. For 
some chemical reaction systems, such as toric dynamical systems (including 
Examples 11.11 and II. 2p . this property is sufficient for proving persistence 
[21 m El [20] , and what was offered in this section are elegant and efficient 
algebraic tools for deriving such proofs. 

Lemma 3.4. For a chemical reaction network G, a siphon Z satisfies prop- 
erty (^) if and only if Z is not relevant (which is equivalent to 55 C "^z)- 

Proof. The "only if" direction is clear. For the "if" direction, let Z he a non- 
relevant siphon. As usual, for a := dimLgtoio we fix a matrix A G 
whose rows span Lcons > and we identify Q with the cone spanned by the 
columns Oj of A. Let F be a facet of Q that contains the image of V>q{^z), 
and let v G M^^^ be a vector such that the linear functional {v, — ) is zero on 
F and is positive on points of Q outside of F. The vector / := vA is in Lcons j 
and we claim that this is a non-negative vector as in ("^). Indeed, li = {v, ai) 
is zero if i G F and is positive if i ^ F, and thus, supp(/) = F'^ Z. □ 

The following result extends Theorem 2 in Angeli et al. [5j. 
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Theorem 3.5. None of the siphons of the network G are relevant if and 
only i/ SatCrc 5S) = (1) if and only if all siphons satisfy property (i^). In 
this case, none of the invariant polyhedra V^(o) has a boundary steady state. 

Proof. The first claim follows from Lemma 13.41 above. The second claim 
follows from the definition of relevant siphons and Lemma |2.2[ □ 

We next present a characterization of the ideals *B and ^^(o) in terms 
of combinatorial commutative algebra. This allows us to compute these 
ideals entirely within a computer algebra system (such as Macaulay 2 [15]), 
without having to make any calls to polyhedral software (such as polymake). 
We assume a subroutine that computes the largest monomial ideal contained 
in a given binomial ideal in the polynomial ring M[ci, . . . , c^]. Let Tgtoi and 
'^cons denote the lattice ideals associated with the subspaces Lgtoi and Lcons- 
These ideals are generated by the binomials c"+ — c"~ where u = n+ — u_ 
runs over all vectors in that lie in the respective subspace. Here, U-^ G Z>q 
and U- € Z2,q denote the positive and negative parts of a vector n in Z'^. 

Algorithm 3.6. The ideals 55 and 5S^{o) can be computed as follows: 

1. The squarefree monomial ideal 53 is the radical of the largest monomial 
ideal contained in Tstoi + (ciC2 • • • c^). 

2. The squarefree monomial ideal 5S^(o) is Alexander dual to the radical of 
the largest monomial ideal contained in the initial ideal ing{o) (Tcons)- 

3. Ifc^^^ is generic (i.e. the polyhedron V^(o) is simple) then the radical of 
in^(o) (Tcons) is a monomial ideal, and its Alexander dual equals 5S^{o). 

The correctness of part 1 rests on the fact that the zero set of the lattice 
ideal Tstoi is precisely the affine toric variety associated with the cone Q. 
Adding the principal ideal (ciC2 ■ ■ ■ c^) to Tstoi is equivalent to taking the 
image of Tstoi in R- The nonnegative variety of the resulting ideal is the 
union of all faces of M>q whose image modulo Lstoi is in the boundary of Q. 

For parts 2 and 3 we are using concepts and results from the textbook 
|18j . The key idea is to use the initial concentration vector c^^^ as a par- 
tial term order. Initial ideals of lattice ideals are discussed in \18\ §7.4]. 
Alexander duality of squarefree monomial ideals is introduced in |18| §5.1]. 
The correctness of part 3 is an immediate corollary to \18\ Theorem 7.33], 
and part 2 is derived from part 3 by a perturbation argument. In the next 
section, we demonstrate how to compute all these ideals in Macaulay 2. 
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4 Computing siphons in practice 



We start with a network that has both relevant and non-relevant siphons. 
This example serves to illustrate the various results in the previous section. 

Example 4.1. We return to the chemical reaction network in Example ll.il 
The sums C+D+E and A+2B+D+2E are both constant along trajectories. 
Chemically, this says that both the total amount of free and bound forms of 
the ligand and the total amount of the free and bound forms of the receptor 
remain constant. Thus, the matrix A can be taken to be 

A / \ / 1 1 1 \ 

A = [aA,aB,ac,aD,aE) = (i2012/' ^' 

The two rows of A form a basis of i^cons- The cone Q is spanned by the 
columns of A. The chamber decomposition of Q is depicted in Figure [H We 
see that the two facets of Q define the following ideal oi Q[A,B,C,D,E]: 

*B = (C, D, E)r\{A, B, D, E) = {AC, BC, D, E) . 

The relevant siphons are derived from 50^^ = {A'^C, AD, E, BC) as follows: 

Sat(imG,5S) = {A,BC,E) = {A,B,E) n {A,C,E) . 

Thus two of three minimal siphons in Example 11.11 are relevant. The third 
siphon is not relevant as its ideal {C,D,E) contains !B. This corresponds to 
the fact, seen in Figure [H that the vectors aA and as span a facet of Q. 

The chamber decomposition of Q consists of three open chambers 0(1), 
0(2), and 0(3), and two rays 0(12) and 0(23) between the three chambers. 
These five chambers are encoded in the following ideals, whose generators 
can be read off from the vertex labels of the polyhedra Vq in Figure [U 

5Sf,(i) = {CD, CE, AC, BC) , 

5Sf,(i2) = {D, CE, AC, BC) , 

»n(2) = (AD, BD, DE, CE, AC, BC) , 

»Q(23) = {AD, BD, E, AC, BC) , 

»f,(3) = {AD, BD, AE, BE, AC, BC) . 

For each chamber O, the ideal SatijMc, ^n) reveals the 0-relevant siphons. 
We find that {A,B,E) is 0(1)- and 0(12)-relevant, and that {A,C,E) is 
0(12)-, 0(2)-, 0(23)-, and 0(3)-relevant. These two siphons define a unique 
vertex steady state on each invariant polyhedron V^io) ■ Note that the vertices 
F^A,B,E} ^-iid F^A,c,E} coincide for polyhedra along the ray 0(12). □ 
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Figure 1: The chamber decomposition of the cone Q for the network in 
Example 11.11 The cone is spanned by the columns of the matrix A in ([5]) . 
Each of the three maximal chambers ^^(1), 0(2), and r2(3) contains a picture 
of the corresponding 3-dimensional polyhedron P^(o) • The vertices of each 
polyhedron are labeled by their supports. The star 'V indicates the unique 
vertex steady state, which arises from the siphon {A, or {^, C, i?}. 
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The need for efficient algorithms for computing minimal siphons has been 
emphasized by Angeli et al. [5], who argued that such an algorithm would 
allow quick verification of the hypotheses of Theorem 13.51 Cordone et al. in- 
troduced one algorithm for computing minimal siphons in We advocate 
Theorem 13. II as a new method for computing all minimal siphons, and Algo- 
rithm [321 as a direct method for identifying relevant siphons. Rather than 
implementing any such algorithm from scratch, it is convenient to harness 
existing tools for monomial and binomial primary decomposition \\2\ I13j . 
We recommend the widely-used computer algebra system Macaulay 2 [15] . 
and the implementations developed by Kahle [16] and Roune [19| . 

In what follows we show some snippets of Macaulay 2 code, and we dis- 
cuss how they are used to compute (relevant) minimal siphons of small net- 
works. Thereafter we examine two larger examples, which illustrate the effi- 
ciency and speed of monomial and binomial primary decomposition. These 
examples support our view that the algebraic methods of Section 3 are com- 
petitive for networks whose size is relevant for research in systems biology. 

Example 4.2. The following Macaulay 2 input uses the command decompose 
to output the minimal primes for the three examples in the Introduction. 

— Example 1 . 1 

Rl = QQ[A,B,C,D,E] ; 

M = ideal(A"2*C, A*D, E, B*C) ; 

decompose (M) 

— Example 1 . 2 

R2 = QQ[e,i,p,q,r,s] ; 

I = ideal (s*e-q, q-p*e, q*i-r) ; 

decompose (I + ideal product gens R2) 

— Example 1 . 3 

R3 = QQ[E,F,P,S_0,X,Y] ; 

J = ideal (E*S_0-X, X*(E*P-X), F*P-Y, Y* (F*S_0-Y) ) ; 
decompose (J + ideal product gens R3) 

By Theorem 13.11 the minimal siphons can be read off from the primes. □ 

Example 4.3. We return to the chemical reaction network of Examples ll.il 
and 14.11 The following Macaulay 2 code utilizes item 2 in Algorithm 13.61 

— Example 1.1: cO-relevant siphons 
cO = {0,0,1,1,0}; 



14 



R = QQ[A,B,C,D,E, Weights => cO] ; 

IG = ideal (A'2*C-A*D, A*D-E, E-B*C, A*B*C*D*E) ; 

ICons = ideal (C*D*E-1, A*B~2*D*E~2-1) ; 

BcO = dual radical monomialldeal leadTerm ICons; 

decompose saturate(IG,BcO) 

In the first line, the vector c*^'^^ was chosen to represent a point in the chamber 
f](l), so the output is the unique r2(l)-relevant minimal siphon. □ 

The next example is of a large strongly-connected chemical reaction, and 
the computation shows the power of monomial primary decomposition. 

Example 4.4. Consider the following strongly connected network which is 
comprised of s species, s — 1 complexes, and s — 2 reversible reactions: 

CiC2 ^ C2C3 ^ C3C4 ^ • • • ^ Cs_iCs . 

The number of minimal siphons satisfies the recursion N{s) = N{s — 2) + 
N{s - 3), where N{2) = 2, iV(3) = 2, and A^(4) = 3. For s = 50 species we 
obtain A^(50) = 1,221,537. The following Macaulay 2 code verifies this: 

s = 50 

R = QQ[c_l. .c_s] ; 

M = monomialldeal apply (1 .. s-1 , i->c_i*c_ (i+1) ) ; 
time betti gens dual M 

We now explain the commands that are used above. First, M is the monomial 
ideal TIq generated by complexes, and dual outputs its Alexander dual [18] . 
which is the monomial ideal whose generators are the products of the species- 
variables in any minimal siphon. Secondly, betti applied to gens dual M 
outputs the degrees of all the generators of dual M; these degrees are exactly 
the sizes of all minimal siphons. The command time allows us to see that 
the computation of the minimal siphons takes only a few seconds. Displayed 
below is a portion the output of the last command above; the list tells the 
number of minimal siphons of each possible size. 



o5 = total 

1 
2 

23 



1 

1 1221537 
1 
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24 



26 
2300 



25 
26 
27 
28 
29 
30 
31 
32 



245157 
497420 
352716 
77520 



42504 



3876 
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The current version of dual in Macaulay 2 uses Roune's implementation of 
his slice algorithm [1^ . For background on the relation of Alexander duality 
and primary decomposition of monomial ideals, see the textbook [18]. □ 

Our final example aims to illustrate the computation of minimal siphons 
for a larger network with multiple strongly connected components. 

Example 4.5. We here consider a chemical reaction network G with s = 25 
species, 16 bidirectional reactions and 32 complexes. The binomials repre- 
senting the 16 reactions are the adjacent 2 x 2-minors of a 5 x 5-matrix (cy), 
and Zg is the ideal generated by these 16 minors CjjCj+ij+i — Cjj+iCi+ij-. 
For this network, Lstoi is the 16-dimensional space consisting of all matrices 
whose row sums and column sums are zero, and Q is a 9-dimensional convex 
polyhedral cone, namely the cone over the product of simplices A4 x A4. 

What follows is an extension of the results for adjacent minors of a 4 x 4- 
matrix in [11, §4]. The ideal 2g is not radical. Using Kahle's software [16] . 
we found that it has 103 minimal primes, of which precisely 26 contribute 
minimal siphons that are relevant. Up to symmetry, these 26 siphons fall 
into four symmetry classes, with representatives given by the following: 

Zl = {Cl4, C21, C22, C23, C24, C32, C34, C42, C43, C44, C45, C52} , 

Z2 = {Cl4, C21, C22, C23, C24, C33, C34, C35, C41, C42, C43, C53} , 

Z'i = {Cl4, C24, C31, C32, C33, C34, C42, C43, C44, C45, C52} , 

Z4, = {Ci4, C24, C31, C32, C33, C34, C43, C44, C45, C53} . 

Under the group of reflections and rotations of the matrix (cjj), the 
orbit of Zl consists of two siphons, and the orbits of Z2, Z^, and Z^ each 
are comprised of eight siphons. The corresponding four types of minimal 
primes have codimensions 13, 12, 12, and 12, and degrees 1, 2, 3, and 6. 

By randomly generating chambers, we found that, for every integer r 
between and 26, other than 23 and 25, there is a point c^''^ in Q such that 
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the number of c*-'^^ -relevant siphons is precisely r. We briefly discuss this for 
three initial conditions. First, let c^^^ be the all-ones matrix. Then V^"^ is the 
Birkhoff polytope which consists of all non-negative 5 x 5-matrices with row 
and column sums equal to five. In this case, all 26 minimal siphons are c^^^- 
relevant: Zi defines a vertex, Z2 and Z3 define edges, and Z4 defines a three- 
dimensional face of "PgC)- Next, consider the following initial conditions: 



V 



1 
1 

1 - e 
1 
1 



1 \ 
1 
1 
1 

1 / 



and 



={0) 



/ 1 
1 

1 

1 

V 1 



1 
1 

1 + e 
1 
1 



1 \ 

1 
1 
1 

1 / 



where e > 0. Again, all 26 minimal siphons are d^'^) -relevant, and Fz^ is a 
vertex, and F^g are edges of V^{o) , but now Fz^ is a five-dimensional 
face. Finally, for initial condition e^'^^ only two minimal siphons are e^°^- 
relevant, both in the class of Zi, and they define vertices. Thus, using the 
results of |3], we can conclude that the system ([1]) is persistent for e^'^). □ 



5 Conclusion 

In the present paper, we gave a method that computes siphons and de- 
termines which of them are relevant. To our knowledge, this is the first 
automatic procedure for checking the relevance of a siphon. As noted by 
Angeli et al. [5], such a procedure is desirable for verifying whether large 
biochemical reaction systems are persistent. Recall that persistence is the 
property that no species concentration tends to zero. In practice, this cor- 
responds to the observed behavior that a substrate that is present at the 
beginning of an experiment will also be present in some amount for all time. 
There is a class of systems for which the non-relevance of all siphons is 
a sufficient condition for such a system to be persistent, so our procedure 
can be used to verify quickly that a large network is persistent. Such a 
class consists of toric dynamical systems [9] . Mathematically, the claim that 
toric dynamical systems are persistent is the content of the Global Attractor 
Conjecture, and we speculate that an algebraic approach to understanding 
siphons may be a step toward the conjecture. 
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